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^^ ' Abstract. We present a simple (mieroscopic) model in which bulk viscosity plays a role in explaining 

the present acceleration of the universe. The effect of bulk viscosity on the Friedmann equations is 
C\| ' to turn the pressure into an "effective" pressure containing the bulk viscosity. For a sufficiently large 

^ I bulk viscosity, the effective pressure becomes negative and could mimic a dark energy equation of 

state. Our microscopic model includes self-interacting spin-zero particles (for which the bulk viscosity 
is known) that are added to the usual energy content of the universe. We study both background 
equations and linear perturbations in this model. We show that a dark energy behavior is obtained 
for reasonable values of the two parameters of the model (i.e. the mass and coupling of the spin-zero 
^— ^ ■ particles) and that linear perturbations are well-behaved. There is no apparent fine tuning involved. 

We also discuss the conditions under which hydrodynamics holds, in particular that the spin-zero 
particles must be in local equilibrium today for viscous effects to be important. 
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1 Introduction 

Experimental evidences — such as type lA supernovae luminosity [1, 2], Cosmic Microwave Back- 
ground (CMB) peaks [3], Baryon Acoustic Oscillations (BAO) [4] and Large Scale Structures (LSS) [5] 
— point toward a universe made of 70% "dark energy" (for reviews see [6, 7]). This dark energy is re- 
sponsible for the observed acceleration of type lA supernovae and is thus repulsive (i.e. it has negative 
pressure). Explaining this acceleration has become a important challenge for theoretical physics. 

There exist many candidate mechanisms/models to explain this acceleration (for reviews see 
[6, 8-10]). The most famous one is to add a cosmological constant to Einstein's equations. Although 
being extremely economical, this solution suffers from severe fine tuning problems. This lead theorists 
to investigate alternative scenarios, such as modified matter content models (e.g. quintessence [11]), 
modified gravity models (e.g. scalar-tensor theories [12-14]) and models relying on the statistical 
distribution of matter in the universe (e.g. backreaction [15-17] or voids [18, 19]). The present status 
is that none of these models offer a satisfactory solution to the dark energy problem because they 
either suffer from severe fine tuning problems, lead to instabilities, arc ruled out or are awaiting 
further analysis. 

In this paper we present an alternative (microscopic) model of dark energy that includes the 
effects of bulk viscosity. In an expanding system, relaxation processes associated with bulk viscosity 
effectively reduce the pressure as compared to the value prescribed by the equation of state. For a 
sufficiently large bulk viscosity, the effective pressure becomes negative and could mimic a dark energy 
behavior. 

The idea of having the bulk viscosity drive the acceleration of the universe is mentioned in 
ref [20]. Unified Dark Matter models take this idea one step further and use a single (bulk) viscous 
fluid to explain both dark matter and dark energy. These studies postulate an exotic equation of 
state for viscous matter of the form C = CoP™ (where Coj"z are parameters) [21, 22] and investigate 
its effect on the evolution of the universe [23-26]. It is interesting to note that the Chaplygin gas 
model [27] is a special case of these bulk viscous models. The conclusion is that it is difficult to 
explain both dark matter and dark energy using a single viscous fiuid [6, 25, 26]. The reason for this 



is that bulk viscosity becomes important at low redshifts (in order to have a negative pressure) and 
density fluctuations in the viscous fluid are quickly damped. The resulting decay of the gravitational 
potential (obeying the Poisson equation) has a dramatic effect on structure formation and leads to 
an integrated Sachs- Wolfe effect on CMB anisotropics that exceeds observational bounds [25, 26]. 

We take a different, more microscopic approach here (for a different microscopic model with 
bulk viscosity induced by dark matter annihilation, see references [28, 29]). Our model consists in 
adding a new component to the energy content of the universe coming from a self-interacting scalar 
field. Thanks to recent progress in thermal field theory, the bulk viscosity of scalar theories has been 
computed from first principles [30]. The resulting functional form depends on physical properties of 
the scalar particle (mass and self-coupling) and is very different from previous viscous models. We 
study the consequences of this new scalar component (including bulk viscosity) on the evolution of the 
universe. We emphasize that our model aims at explaining dark energy only and the arguments against 
Unified Dark Matter models with bulk viscosity do not apply here. Indeed, the scalar component of 
our model is added on top of the usual Cold Dark Matter (CDM) component and thus the late 
damping of fluctuations in the viscous fluid just results in dark energy being very homogeneous today 
without affecting structure formation or the CMB. 

The rest of the paper is organized as follows. In section 2 we present some theoretical concepts 
related to hydrodynamics and bulk viscosity. We also present our model and the conditions under 
which it is valid. Sections 3 and 4 are dedicated to the study of background equations and linear 
perturbations respectively. In each case we derive the modifications to the equations due to bulk 
viscosity and show the results of simulations done with the Cosmic Linear Anisotropy Solving System -'^ 
(CLASS) code [31, 32]. The compatibility with observations of our model is presented in section 5. We 
then conclude in Section 6. 

2 Theoretical concepts and presentation of the model 

Hydrodynamics is an effective theory that describes the evolution of "fiuid cells" , where the equations 
of motion are given by exact local conservation of energy-momentum (we implicitly assume that there 
are no other conserved charges in the following). A fluid cell is a macroscopic average (of size I) over 
many fundamental fluid constituents. For a perfect fluid, macroscopic quantities such as pressure 
p, energy density p and velocity v are smooth inside a fluid cell and local thermal equilibrium is 
maintained. For a viscous fluid, microscopic interactions make p, p and v vary appreciably over a 
mean free path and local thermal equilibrium is not maintained inside a fluid cell. The effect of 
dissipation due to microscopic interactions is characterized by transport coefficients (shear and bulk 
viscosities). In the following we focus on bulk viscosity since we are interested in the universe and 
there is no shear in an expanding system. 

Viscous corrections to the homogeneous evolution of the universe are discussed in great details 
in ref. [33]. Let us repeat the salient points. The energy- momentum tensor for a perfect fluid is given 
by: 

T^:,,,,,^P9^'' + {p + pWU\ (2.1) 

where g^^^ is the metric tensor and t/^ is the velocity four- vector (normalized such that U^U^ = — 1). 
In a Friedmann-Robcrtson- Walker universe with metric g^'^ — diag(— l,a, a,a), the usual procedure 
leads to the Fricdmann equation: 



and the energy conservation equation: 
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where a is the scale factor, k is the spatial curvature and dots represent derivatives with respect 
to time. To include viscous corrections, we need the expression for the viscous energy-momentum 
tensor. Transport coefScients characterize linear deviations away from equilibrium due to thermody- 
namic forces (gradients). In a viscous fluid, momentum flows between neighboring fluid cells due to 
microscopic interactions and is thus driven by velocity gradients. The most general energy-momentum 
tensor linear in velocity gradients that satisfies the second law of thermodynamics for all fluid config- 
urations is [33]: 

r;4ou. = TpTrfcct - C ig'- + U^Un D,U\ (2.4) 

where C is the bulk viscosity and D-^ is the covariant derivative. In a Friedmann-Robertson- Walker 
universe, we have D-yV^ = Za/a and we re-write the viscous energy-momentum tensor as: 

^vILu. = pu'u'^ +(p- 3C^) {g^'- + u^un 

= P,sg'"' + {p.ii + p)U^U\ (2.5) 

From eqs 2.1 and 2.5, we see that the effect of bulk viscosity is to change the pressure p to an 
effective pressure Peff = P — 'ii^a/a. The physical interpretation is clear. An expanding fluid leaves 
its equilibrium state; the energy density decreases and the pressure also decreases. In the absence of 
bulk viscosity, the fluid relaxes instantaneously and pressure and density are related by the equation 
of state. Bulk viscosity dampens this behavior by introducing a finite relaxation timescale, hence 
producing a shift between the equation of state pressure and the true pressure. 

The background evolution equations are obtained by replacing p by Poff in the energy conservation 
equation 2.3 (the Friedmann equation 2.2 remains unchanged). We note that for a large enough C,, 
the effective pressure becomes negative and could mimic a dark energy behavior. The bulk viscosity 
coefficient C, depends on microscopic interactions and must be computed from a more fundamental 
theory. 

2.1 Presentation of the model 

The background equations 2.2 and 2.3 are effective equations that do not rely on microscopic physics, 
except through the bulk viscosity parameter. We thus need some microscopic physics to "personify" 
the bulk viscosity. For reasons we discuss in section 2.3, no Standard Model particles can produce a 
bulk viscosity in the late universe. We thus need to add a new component to the energy content of 
the universe in order to induce some bulk viscosity in the present universe. 

For simplicity, we add a massive scalar field with quartic interactions to the Lagrangian of the 
Standard Model (cubic interactions would change our results in a trivial way and are left out in the 
following) : 

L - \d,(pd>^<t> - Imlcf" - ^<p\ (2.6) 

We also make several assumptions. First we require that the assumptions of hydrodynamics for this 
scalar fluid are satisfled. The validity of hydrodynamics is crucial for the concept of bulk viscosity to 
make sense. These assumptions include local thermal equilibrium, small inhomogcneities (gradients) 
and absence of instabilities. Each of these points arc discussed in section 2.3. We also assume that 
the scalar and Standard Model particles are in thermal equilibrium in the early universe and decouple 
from each other afterward. After this decoupling time, the temperature of the self-interacting scalars 
has a distinct evolution from that of photons and other Standard Model particles in equilibrium 
with the photons. Between electron-positron annihilation and today, the photon temperature scales 
as T-y oc a^^ . We shall see later that the scalar temperature T, scales in the same way at least 
between the time of scalar-photon decoupling and the non-relativistic transition (Ts ~ toq). Hence, 
the ratio T^jT^ is expected to remain constant over a long time interval, typically from Big Bang 
Nucleosynthesis until the end of radiation domination. 



The model has three free parameters: the mass ttiq, the self-couphng A of the scalar particles, and 
the ratio of scalar to photon temperatures e = Tg/T^ after electron-positron annihilation. Note that 
e < 1 because (aTj) can increase due entropy releases until electron-positron annihilation, while (aTs) 
stays constant (no interaction with the Standard Model particles) after scalar-photon decoupling. The 
limit e = 1 correspond to late decoupling not followed by any entropy release and the limit e <C 1 
corresponds to early decoupling with significant entropy releases. 

We stress that our model is very different from quintessence models. Both include a scalar field 
with a potential, but the range of parameters is different. In quintessence models, the scalar field 
is seen as classical and coherent over large scales, and the potential is adjusted in such a way as to 
obtain slow-roll conditions in the late universe. In contrast, we are dealing here with a fluid of scalar 
particles, and we adjust the parameters of our model in order to allow collective (bulk viscous) effects 
to play a significant role in this imperfect fluid in the present universe. 

We also mention that our model is not a generalization of previous bulk viscous or Chaplygin 
gas models. The functional form we use for the bulk viscosity (see section 2.2) is very different from 
the (generalized) Chaplygin form. The bulk viscosity for a scalar theory has been computed from first 
principles and depends on the physical properties (mass and self-coupling) of the scalar particles; no 
such interpretation is available for the (generalized) Chaplygin gas form. 

2.2 Pressure, density and bulk viscosity of a scalar field 

We need the pressure, energy density and bulk viscosity of a scalar field to integrate eqs. 2.2 and 2.3. 
The pressure of a gas of massive scalar particles (with vanishing chemical potential) is given by (e.g. 

[34]): 

P=^^^{^)+OiXTfl (2.7) 

where K2 is a modified besscl function of the second kind. In the following we neglect 0{\T^) since 
we expect bulk viscous effect to dominate at low temperature. Using the thermodynamic relation 
p = T^(? {p/T) / dT, we obtain the corresponding energy density: 
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In a remarkable tour de force, Jeon has computed the leading order bulk viscosity in a scalar 
theory from diagrammatic methods [30]. The computation involves the resummation of an infinite 
number of diagrams and can only be obtained numerically. A best fit of the numerical result is: 

where the constants ki, K2 and K3 are 12.13, 1.54 and 1.13, respectively. The thermal mass (i.e. 
the sum of the zero-temperature mass plus thermal modifications) is given by rn^Yi — "*o + '^"^thi 
where the thermal modification is Sm^y^ = XT^/24. The quantity rh^ = m"^^ — [Tg /2){dm1^/ dTg) 
gives the breaking of scale invariance in the theory (at lowest order). More precisely, it is given by 
m"^ = ml- l3{X)T^/48 where the beta function is /3(A) = 3A2/(167r2). Note that eq. 2.9 is only vahd 
for weak coupling A < 1. 

We can understand eq. 2.9 in the following way. Bulk viscosity characterizes the relaxation of an 
expanding fiuid. As a result, it is proportional to the mean free time between collisions of the fluid 
constituents. The mean free time is itself inversely proportional to the cross section of the processes 
responsible for re-establishing equilibrium. In an expanding system, equilibrium is achieved by number 
changing processes (2 — >■ 4 in Xcj)"^). This explains the A~^ dependence of C. The exponential is due 
to Boltzmann suppression of final states for these number changing processes. 

The bulk viscosity also vanishes in a conformal theory, since in that case dilatation is a symmetry 
and the fluid cannot leave equilibrium. It is thus proportional to the (square) of the breaking of scale 
invariance (i.e. to^). Consequently, C. has a very different coupling constant behavior at different 
temperatures. At low temperature m'* ~ ttiq and ^ ^ A^^; at high temperature m^ ~ (3'^ and 
(^ ^ Xy'^^~^>/'^, implying that the bulk viscosity is highly suppressed in this regime. 



2.3 Validity of the hydrodynamic approximation 

The crucial point in this model is that the new scalar degree of freedom can be considered as a 
fluid in the dark energy domination epoch. Otherwise, the concept of bulk viscosity does not make 
sense. Thus the model must satisfy the assumptions of hydrodynamics. Let's study each assumption 
separately. 

Local equilibrium. For a hydrodynamic description to be valid, the fluid constituents must inter- 
act sufficiently often so as to maintain (approximate) local equilibrium (the constituents must 
not be "frozen out"). More technically, the need for local equilibrium is reflected in the fact 
that transport coefficients are linear deviations away from equilibrium and are computed from 
equilibrium correlators [30]. 

The criterion to be in equilibrium in an expanding universe is Tscaim/H 3> 1, where Fgcaiar is the 
rate of processes that are ultimately responsible for the bulk viscosity and H is the expansion rate 
of the universe. The rate Fgcaiar is estimated as follows [30]. First we have that Fgcaiar ~ l/Tfree 
where Tfrcc = n/{dW/dVdt) is the mean free time between collisions, n is the density of particles 
and dW/dVdt is the transition rate per unit volume. The processes that are responsible for bulk 
viscosity in A;/)^ theory are 2 — >■ 4 scatterings. We expect the bulk viscosity to dominate at low 
temperature due to the exponential factor in eq. 2.9. For low temperatures T <^ rriQ we have 
n = dN/dV - O(e-™o/^-) (Boltzmann suppression) and dW/dVdt ~ 0{X'^e-^"'°/^=) (there 
should be a rriQ to get the units right). Thus the rate of bulk viscosity processes are: 

r.ealar ^ - ^ n^oX^e-^o/T. ^ ^^^,^-3n.o/T. ^ (2.10) 
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The equilibrium criterion becomes: 



H 



> 1. (2.11) 



This criterion must be fulfilled until the present epoch for the model to make sense; this is thus a 
constraint on the parameters of the model. This also explains why no Standard Model particles 
produce a large bulk viscosity in the late universe, since they all decouple before dark energy 
domination. 

Small inhomogeneities and gravitational instabilities. Hydrodynamics is an expansion in gra- 
dients (see for example eq. 2.4). For this expansion to be well defined, velocity gradients must 
be small (small inhomogeneities) . Similarly, clumps of scalar matter may form due to gravity. 
If these clumps become too large (exponential instabilities) , then the fluid breaks apart and a 
hydrodynamic description is not appropriate anymore. 

It is difficult to estimate the size of gradients a priori. For that, a full numerical study of linear 
perturbations is necessary (a dynamical system analysis might also shed some light on these 
issues, e.g. [35]). We perform such a numerical analysis in section 4. Our results show that 
linear perturbations are well behaved and no instability develops (see section 4.2). 

It is important to emphasize that we are using first order (in gradients) relativistic hydrodynam- 
ics. It has been pointed out in Ref. [36] that first order hydrodynamics is acausal and unstable 
(for certain modes); to cure these problems, second order hydrodynamics should be used instead 
(see [37] for a cosmological case study). A more recent analysis shows that acausality implies 
unstable behavior [38, 39]. The physical explanation is that an acausal disturbance propagating 
in a light cone would create a singular behavior at the edge of the light cone, since a light cone 
cannot be crossed in a covariant theory. In our case, growing modes are damped out by bulk 
viscosity and no instability develops. We conclude that our hydrodynamic equations are well 
defined and do not produce acausal signal propagation. It thus seems that the use of linear 
relativistic hydrodynamics is justified for the modes considered here. 



Thermodynamical instabilities. This instability is called "cavitation" (for a discussion, see e.g. 
[40]). In everyday engineering systems, cavitation manifests itself when the pressure of a fluid 
goes below its saturated vapor pressure. In such a case, the fluid becomes unstable to the 
formation of bubbles. The formation of bubbles is akin to the formation of clumps; if there are 
too many bubbles in the fluid, then a hydrodynamic description is not suitable anymore. 

The process of bubble formation is called nucleation. There are two types of nucleation: homo- 
geneous and heterogeneous. Homogeneous nucleation happens when a fluid is pure and smooth. 
Heterogeneous nucleation happens when a fluid has impurities and is in contact with boundaries. 

In our model, the saturated vapor pressure is the vacuum pressure (i.e. zero). Thus when 
the pressure becomes negative (required for dark energy domination), the fluid becomes in 
principle unstable to the formation of vacuum bubbles. In the following we use the language of 
homogeneous nucleation since the fluid is very pure (except maybe for "gravitational" impurities) 
and smooth (no boundaries). 

In everyday engineering systems, it is possible to reduce the pressure of a fluid under its saturated 
vapor value without forming bubbles (a bit like a supercooled fluid). In other words, a normal 
fluid can withstand "some" under saturated vapor pressure, it is not automatic that bubbles 
form. The physics behind that behavior is the following. Under normal conditions, there is 
formation of voids between the fluid constituents due to thermal fluctuations. When a fluid 
is subjected to a pressure that is under the saturated vapor pressure, these voids grow and 
form bubbles. But there is a counter force to this growth due to the potential between fluid 
constituents. When the outward pressure is greater than the force between constituents, bubbles 
grow; if not, then the fluid stays as it is. This property of a fluid is called "tensile strength". 
For homogeneous nucleation, tensile strength depends only on the properties of the fluid; in 
heterogeneous nucleation, this tensile strength is strongly dependent on impurities and boundary 
conditions. 

Let's denote pt max the maximum pressure our scalar fluid can sustain without breaking apart 
(i.e. its tensile strength) and puE the negative pressure that makes the tiny voids due to thermal 
fluctuations grow. In the dark energy domination epoch, this pressure is given by pde = wpDE, 
where w is the effective equation of state and pde ^ 10~^^ eV* is the present dark energy 
density. A criterion for the absence of cavitation can be phrased in the following way: there is 
no cavitation in the system if the tensile strength is larger than the maximal outward pressure: 

Pbmax > I^cIPDE- (2.12) 

This condition must be satisfied for a sufficiently long time in order to have dark energy dom- 
ination without the fluid to break apart. The tensile strength ph max is a property of the fluid 
and may depend on A, m and the dimensionless inverse temperature x = toq/Ts- Thus the 
above criterion for no cavitation is a constraint on the parameters of the model. 

We can roughly estimate the scalar fluid tensile strength using the theory of homogeneous 
nucleation (e.g. [40]). The result is that the tensile strength is given by: 



where S is the surface tension and J is the nucleation rate. Since the above relation depends 
logarithmically on the nucleation rate, we expect a very weak dependence on J. Thus for 
estimate purposes, we set the logarithm to one. The surface tension is a force per unit length; 
we thus expect it to be proportional to the coupling A and a quantity that has dimensions of 
[Energy]"^. The two relevant energy scales in the problem are tuq and T^. From dimensional 
analysis, we obtain four possible forms for the surface tension: S ^ ATj*, S ^ XuioTg, S ~ Atti^Ts 
and S ~ Attiq. The first one is expected to be valid only at ultra-relativistic energies and we 
discard it. The last one is unrealistic because of its lack of temperature dependence. Empirical 



laws show that the typical behavior of surface tension is linear in temperature. We thus use the 
form S ~ Xm^Ts in the following analysis. Plugging S ^ Xm^Ts in Eq. (2.13) we get: 

.._-(—) ^r^. (2.14) 

Combining this with the no cavitation criterion, we obtain the following constraint on the model 
parameters: 

Cavitation is also studied in the context of rclativistic heavy ion collisions [41]. In this study, it 
is stated that any hydrodynamic simulations should be stopped when the pressure goes below 
zero because of cavitation. We argue here that, just like in other engineering systems, a fluid 
can sustain some negative pressure without cavitating. 

3 Background evolution 

3.1 Implementation in CLASS 

The pressure, energy density and bulk viscosity arc all functions of the dimensionless inverse temper- 
ature X = mo/Ts (c.f. eqs. 2.7-2.9): 

P.W - ^^. (3.1) 
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where primes denote derivatives with respect to x and the subscript s relates to the scalar particles. 
The energy conservation equation relates ps to ps and Pcs = Ps — 37?^. The ususal procedure is 
to use the equation of state to express the pressure in terms of p and then solve for p. In our case, the 
results for the pressure, energy density and bulk viscosity arc all expressed in terms of x. To integrate 
the equations, we re-write the LHS of the energy conservation equation 2.3 as {dps/dx){dx/dt), where 
the derivative dps/dx is known. The resulting background equations requiring a numerical integration 
over time are: 

a = aH , (3.4) 

where H is inferred from the Friedmann equation, including a contribution p^ to the total energy 
density, and p'^ is given by: 

p;(x) = ?4f-^ + ^-^^V (3.6) 

27r^ \ x-' x'^ X J 

Equations 3.4 and 3.5 are integrated using the Cosmic Linear Anisotropy Solving System (CLASS) 
code [31, 32] modified to take into account the new scalar energy component. 

The minimal extension of the flat ACDM model in which the cosmological constant is replaced by 
a viscous scalar fluid has eight free parameters, namely the three fundamental parameters of the scalar 
sector (ttiq, A, e) and five parameters describing the other sectors. These five parameters are chosen 
to be (rif,, ft,, As^ Hs, t) (baryon density fraction, reduced Hubble parameter, primordial spectrum 



amplitude and tilt, reionization optical depth). In principle, each combination of {mo, A, e, h) leads 
to a unique value of the current scalar density fraction fig, from which one can infer the CDM density 
fraction flcdm = 1 — 51^ — Qs- The fact that the universe is flat imposes a bound in the {mo, A, e, h) 
space ensuring that il^ < 1. 

In order to compare our model with observations, it is much more convenient to use fig or 
ri„i = 1 — fi^ as one of the model parameters. In that case, mo can be seen as a function of (fim, A, 
e, h). This paramctrisation is the one that we implement in CLASS. For each model parameters 

(fi,„, rtb, h, As, Us, T, A, e) , (3.7) 

the code searches for the value of mo leading to the correct relic density with a simple bisection 
algorithm. This involves running the CLASS module background. c several times with different ryip's 
until the correct mass is found up to a given accuracy. The other modules are run after the obtention 
of the correct mass. 

3.2 Results for the background evolution 

Figure 1 shows the evolution of the densities, temperature and effective pressure in the minimal model 
including a viscous fluid described by the eight parameters listed above (note that only five of them 
are relevant for the background evolution). Solid lines correspond to the choice 

(^n, ^b, h, A, e) = (0.28,0.05 ,0.72,0.25,0.7), (3.8) 

giving TTiQ = 1.0 eV with a viscous fluid temperature Ts = 0.7 T-y of the same order as the neutrino 
temperature as long as the fluid behaves like radiation. Dashed lines correspond to the same pa- 
rameters as above except for e = 0.1, implying that a number g^ ^ 100 of relativistic degrees of 
freedom annihilate after the scalar-photon decoupling, and transfer their entropy to the photons and 
particles in thermal equilibrium with them. Finally, dotted lines correspond to a model with the same 
parameters as above except for A = 10~^, giving toq = 2.7 x 10~^ eV. 

On the left panel, we sec that scalar species are initially relativistic with ps oc a~^ and Tg oc a^^. 
This remains true as long as a; ^ 1 (Tg ^ mo) and bulk viscosity is negligible. The contribution of 
the scalar fluid to the commonly used "effective neutrino number" Ngg is given by 

AiVcff = i X ^ (^^y '^4 = 2.2e4, (3.9) 

where the factor is obtained from the neutrino-to-photon temperature ratio assumed in the usual 
deflnition of iVcff and from the fact that we are considering bosons (AA'off — 1 corresponds to the 
case of one fermion plus its anti-particle). 

When a:; ~ 1 (T^ ~ m,o), the bosons become non-relativistic and their density starts to dilute like 
that of ordinary matter {ps oc a~^). This is clearly visible for the model with solid or dashed lines in 
figure 1. In the case of coUisionless species like massive neutrinos, the temperature keeps decreasing 
like a~^ during the non-relativistic regime. This is not true for self-interacting species like the scalar 
fluid considered here. An analytic study of the evolution equations in the limit x 3> 1, C = shows 
that the temperature evolves according to Ts e~"^°^'^^ oc a~^ during this stage. 

We also see that the viscosity ( increases with time. At some point, the terms 3H( and Ps 
become of the same order. After that moment, the effective pressure pcs ~ Ps — ^H( changes sign 
and flnally approaches —ps- The system then tends towards a fixed point with constant values of ps, 
Ps, H, (^ and Ts such that Peff = —Ps (and Ps <C 3i/C)- As expected, this mechanism mimics dark 
energy at least at the level of the background evolution. 

3.3 Parameter dependence 

It is interesting to compare the three models displayed in figure 1. In the model with a smaller e 
(dashed lines), the fluid density and temperature are initially smaller (in this case, N^s — 3.04 like 
in the minimal ACDM scenario). Since the mass is unchanged with respect to the reference model 



with solid lines, the transition to the non-relativistic regime takes place at the same temperature 
Ts ~ 1 eV, corresponding to an earlier time. As a consequence of the smaller temperature, viscosity 
becomes important earlier and the scalar fluid starts its dark energy behavior earlier. In the model 
with a smaller value of A (dotted lines), and also a smaller value of toq corresponding to the same 
relic density ils, the early and late behaviors are similar to those of the reference model. However, 
we do not observe an intermediate non-relativistic stage with ps oc a~^, because the time at which 
Ts ^ Too coincides with that at which viscosity becomes important: the fluid goes directly from a 
radiation-like to a dark-cncrgy-likc behavior. 

The relation between A and toq for a given (f7s, h) or (f7„i, h) is difficult to estimate analytically. 
In principle, this relation can be obtained by noticing that close to the fixed point one has: 

p° = f],p° ~ -peff . (3.10) 

The effective pressure is then dominated by the bulk velocity contribution, so we can rewrite these 
two equalities (where the second one is approximate) as: 

p0^n,^c3HoC. (3.11) 

Replacing ps and ^ by their expression in terms of x and other fundamental parameters (c.f. eqs. 
3.2, 3.3), we obtain two independent relations between (a;, toq, A,i7o7 ^s)- These two equations can 
be combined in order to extract x and find the temperature at the fixed point as a function of 
(to-o, a, Hq, fls). Finally, x can be replaced in any of the two equations in order to obtain a relation 
between (toq. A, Hq, fJs). Note that e plays no role in this calculation. 

This system of equation is difficult to solve analytically. However we obtain such a relation 
numerically, thanks to the bisection algorithm implemented in CLASS. Namely, we vary A in the range 
[10"^^°, 10^], h in the range [0.65, 0.75], J7,„ = 1 — ils in the range [0.2, 0.4] and obtain the corresponding 
niQ in each case. We show in the left panel of figure 2 the corresponding region in (A, toq) space (thin 
blue band). This region is compatible with hydrodynamical constraints only in the range A > 2 x 10^^ 
(otherwise the fluid is broken apart by cavitation), corresponding to mo > 0.01 eV. The fluid is always 
in local equilibrium in this range of coupling. The scaling between A and ttiq changes around the value 
A ^ 10"^, but in the physically interesting range A ^ 10~^ we flnd a fltting formula accurate at the 
percent level: 

/ r-, \ 0.295 / . \ 0.445 

„„..«SA»-(i^) (_i_) .V. (3.12, 

Note that the validity of the model for a large self-coupling constant is far from obvious, since the 
expression for bulk viscosity is obtained using perturbation theory and assumes A < 1 . If we impose 
a theoretical prior A < 1, the allowed range for the mass reduces to ttiq G [0.01, 2] eV. 

We checked numerically that the above relation between ttiq and A is not affected by e, as expected 
from our previous reasoning. Moreover, the current ratio x'^ = mo/T^ depends on (A, Hq, fls) but not 
on e. Hence the e parameter is only important at early time through its impact on the total radiation 
density. 

We see from flgurc 1 that the viscous fluid does not reach the flxed point (or stabilize with a 
constant energy density) immediately. This suggests that our model predicts a value of the equation- 
of-state parameter if = Ps/Ps larger than -1 today and slightly varying with time. A priori, the 
current value wq could depend on (A, Hq, fig) and even slightly on e. A careful numerical investigation 
illustrated by the right panel of flgure 2 shows that wq depends only on A and Qs (or flm)- The same 
appears to be true for its time-derivative. This seems to imply that the present model is an example 
of a bulk viscosity- induced dark energy model that does not "cross the phantom divide" [42-44] . 

In summary of this subsection, we flnd that the model under consideration has the same back- 
ground evolution as a minimal ACDM model extended in the following way: 

• It has extra relativistic degrees of freedom at early times, controlled by the parameter e. The 
equivalent effective neutrino number is 

TVoff = 3.04 + 2.26^ . (3.13) 




Figure 1. Left: Evolution of the radiation, matter and viscous fluid densities in a model with h = 0.72 and 
Qm = 0.28. Solid lines correspond to the parameter choice (e, A) = (0.7,0.25) leading to mo — f.O eV. The 
dotted (resp. dashed) line is a variant with (e. A) = (0.1,0.25) (resp. (0.7, 10"*)) leading to mo = 1.0 eV 
(resp. mo = 2.7 x lO"'^ eV). Right: For the same models, evolution of the viscous fluid temperature and of 
the correction to the effective pressure induced by viscosity. 





0.0001 0.001 0.01 



Figure 2. Left: Relation between mo and A in a flat universe with 0.2 < Sim < 0.4 (i.e. with 0.6 < flgoo < 0.8), 
assuming 0.01 < e < 1 and 0.65 < h < 0.75. The shaded region is excluded by constraints from hydrodynamics 
(cavitation). Right: For the same range of parameters e and h, current value of the equation of state parameter 
of the viscous fluid as a function of A, for three different values of Qrn- 

If scalar particles decouple from the rest of the plasma very late (but before electron-positron 
annihilation), e may be in the range [0.5,0.7] and produce a significant increase of Ngff ^. If they 
decouple early enough such that e < 0.5, the increase in iVgff is negligible given the accuracy of 
future CMB and LSS experiments. 

• It incorporates slowly decaying dark energy, so our model mimics the so-called wCDM model [45, 
46] rather than ACDM. The value of w today is found to be controlled by {ils, A) or equivalently 
by {ilrrn A), and the follovifing fit is accurate to better than one percent: 



Wo 



-0.9085 -I- 0.21(0,™ -0.3) 



3^-0.4^0-4 



(3.14) 



^Indeed, if the bosons decouple before electron-positron annihilation, their temperature is at most equal to that 
of neutrinos, T^ = (4/ll)^/''T-y ~ 0.7T-,. In principle, it would be possible to assume that our bosons decouple 
from the plasma after electron-positron annihilation, corresponding to e in the range [0.7,1]. The physics of Big Bang 
Nucleosynthesis and photon recombination would then be affected in a more non-trivial way than through N^g . We do 
not consider this case for simplicity. 
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• For a very precise description, one should take into account the fact that the parameter w 
varies slowly with time during dark energy domination. If we try to capture this small redshift 
variation with a standard first-order expansion of the type w{a) = wq + Wa(l — a/oo) [45, 46], 
we find numerically that the fit 

Wa = 0.129 + 0.02(a„ - 0.3) + A""•^10"^ (3.15) 

is accurate up to 5%. 

It remains to be checked that the equivalence between our model and such an extended wCDM model 
is valid also at the level of perturbations. The study of linear perturbations is the object of section 4. 
Since for A 3> 10"'* the scalar fluid undergoes a non-relativistic regime with ps ex a~^ prior to 
the dark energy regime, one could think that it could play alternatively the role of dark matter and 
dark energy. At the level of the background evolution, this could be achieved by setting Qcdm = 0, 
flm = ris + ^b = 1, and by increasing A (and hence mo) until the scalar density during the non- 
relativistic regime matches the usual CDM density. We find that this is indeed possible, but for 
values of A much larger than one. We conclude that our bulk viscosity model cannot explain dark 
matter and dark energy at the same time, but in principle one could try to extend this model with 
such a goal. 

4 Linear perturbations 

As pointed out in section 2.3, it is important to study linear perturbations around the background 
solution to discuss the issue of small inhomogeneities and instabilities. In the following we derive the 
equations for linear perturbations (including viscous effects) in the synchronous gauge, since it is the 
gauge that is normally used in numerical simulations. 

4.1 Equations for the perturbations 

Our starting point arc the standard evolution equations for scalar perturbations in the synchronous 
gauge (wc follow the conventions of Weinberg [47] ) : 

Sp + 3-i5p + 5p)-KTip + p)6U+ip + p)i'-k'^-7:'^ = 0, (4.1) 

a a'' a 

{p + p)SU+ f{p + p) + 3-{p + pyj6U + 6p-pTr^ = 0, (4.2) 

4> + 2-iP + AirGiSp + Up - Pn'^) = 0, (4.3) 

a 

where p and p are the background values for the density and pressure, Sp and Sp arc the pertubations 
in density and pressure, SUi = diSU + 5UY are the velocity perturbations (where diSUY = 0), ?/! is 
a linear combination of metric perturbations, tt is the anisotropic stress tensor, fc is a wavevector 
and dots represent derivatives with respect to proper time. Equations 4.1 and 4.2 correspond to 
energy conservation and momentum conservation, respectively. There is one set of energy-momentum 
conservation equations for each type of matter that does not exchange energy and momentum with 
other constituents. In the following we present the treatment of linear perturbations for the added 
scalar matter, but it must be borne in mind that we consider all linear perturbations in our final 
analysis. 

By assumption the scalar fluid in our model only interacts gravitationally with other particles. 
Thus we can write a set of energy-momentum conservation equations for the scalar perturbations 5ps 
and 5Us- Since equations 4.1- 4.3 do not depend on the particular form of the energy-momentum 
tensor, we make the replacements p — > pcff and 5p — > 5pcs to account for bulk viscous effects. We also 
assume that the anisotropic stress tensor tt^ for our scalar fluid is zero. There is no data to constrain 
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this number and we make this assumption to simphfy the equations. The evolution equations of 
perturbations for our scalar fluid read: 



5ps + i-iSps + 5pcS) :^{Ps +Pcfi)5Us + [ps +Poff)V' = 0, 

a 



{ps +Pcs)5Us + I [ps +Poff) + 3-(p,, +pcff) 1 5Us + 5pc« = 0, 

ijl + 2-V' + AllG{5ps + (5pothcr + iSpcS + 3(5pothcr) = 0, 



(4.4) 
(4.5) 
(4.6) 



where the subscript "s" refer to the scalar fluid and "other" to the other fluids present (photons, 
neutrinos, baryons, dark matter). We note that all fluid perturbations are necessary to obtain the 
correct metric perturbations. 

In our model, the forms of ps and poff are given in cqs. 3.1- 3.3. Wc note that the energy density 
is only a function of x. So we can express a perturbation in energy density solely as a perturbation 
in X (or cquivalcntly a change in temperature): 



5ps = D{x) 5x, 



(4.7) 



where the function D{x) = p's{x) is given by eq. 3.6. The change in effective pressure is a bit more 
tricky, since Peff also depends on D^U^ . We write: 

(5peff = Sips-CD^V) 

= Sps - 5(:{D^U'<) -CSiD^V) 

= P,s{x)5x^C5{D^U^), (4.8) 

where the function Poff(a;) can be obtained from cqs. 3.1 and 3.3. To obtain 5{D^U'^), we write: 



6{D,U'') = 5 



dx" 
d{5ir)_ 
dx" 



KxU' 



+ r:,{6u') + {6r:,)u' 



^^ + 3-isuO)Hsr:.)u\ 



(4.9) 



where in the last line we used T];,^{6U^) = (r[!^ + Ti^)SU^ = ^qSU" + Tl^SW = {3a/a)6U" when 
spatial curvature is zero. The 4-velocity is [/^ = (— l,w') and thus SU'^ = 0. Keeping only linear 
terms in velocity and perturbations and using the perturbed Christoffel symbols in the synchronous 
gauge [47], we obtain: 



SiD^U") = a-^V'^SU + Tp. 



(4.10) 



The perturbation of the effective pressure can thus be expressed in terms of temperature, velocity 
and metric perturbations: 

Spcs = Pcsix)Sx - C{a-^V^SU + i>). (4.11) 

Using eqs. 4.7 and 4.11, we can re-write the evolution equations for the perturbations 4.4- 4.6 as: 



6x + 



D + 3H{D + Peff) 



D 



Sx - 



3CH 



^jPs+Pcff) 



D 



6Us 



{p,+p,s-3HC) 



{ps +Pcs)SUs + PcsSx - 



D 



^2 1 



^ = 0, (4.12) 



(Ps + Pcff ) + 3H{ps + Pcs) + C^ 



SUs -Ci^^O, (4.13) 
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ip + [4TrG{D + 3Poff )] 5x - 



12ttGC^ 



SUs + [2{H - GnGO] ^j + 4^G((5pothcr + 3(5pothcr) = 0. (4.14) 



These are the evolution equations for scalar perturbations for the viscous (scalar) fluid. We integrate 
them using CLASS and the results are shown in section 4.2. For completness and in order to make 
our results easily reproducible, we provide in Appendix A a complete list of the equations added to 
the CLASS perturbation module. Essentially, these are equations 4.4, 4.5 and 4.11 translated into the 
notations of Ma & Bertschinger [48], which are used throughout the CLASS code. 

4.2 Results for the perturbation evolution 

We first study the evolution of a few characteristic scales related to perturbations in the viscous 
fluids. For a perfect fluid, the qualitative behavior of perturbations is essentially captured by the 
Jeans wavenumber: 

where Ci is the sound speed in the fluid. Note that the sound speed is usually deflned as the variation 
of the pressure with respect to the density for constant entropy. Since a viscous fluid contains entropy 
perturbations, we cannot readily obtain c^ from the ratio Spes/Sps- We should instead bear in mind 
that the viscosity of the fluid remains negligible in an adiabatic transformation, so we must deflne the 
sound speed for the scalar fluid using the equilibrium pressure: 



2 dpcs 

c. - 



dps 



Ir^ii^ 



The Jeans wavenumber obtained with Ci — Cs is shown as a solid blue line in figure 3 (for the 
same reference model as in the previous section). For comparison, the dashed line is obtained using 
c^ = Pcs{x)/D{x), in order to give an indication of when the viscosity becomes important and which 
scales it can affect. 

Let us now summarize the expected behavior of perturbations in the scalar fluid when their 
wavenumbers cross the four different regions shown in flgure 3 (left). After crossing the Hubble radius 
k ~ aH, the modes experience gravitational collapse as long as fc < fcj. This stage is very brief: even 
if the scalar particles become non-relativistic before the dark-energy-like stage (as it is the case in 
our reference model), the sound speed decreases very slowly below l/v3 and the Jeans wavenumber 
remains of the same order of magnitude as aH. When the modes cross the Jeans length (fc > fcj), 
and as long as viscosity is negligible, acoustic waves start to propagate in the fluid with a velocity Cs ■ 
In Fourier space, 6s oscillates with a pulsation kcs- 

At some point, viscosity becomes important, and due to non-adiabatic contributions |(5poff| be- 
comes very large with respect to Cs|<^Ps|- In our reference model, this happens around r = 3000 Mpc, 
since this is the time at which the two characteristic scales plotted in flgure 3 become different. The 
very large value of \6pcff/5ps\ implies that the friction term in the equation of propagation oi Sg drives 
the perturbation to zero, on a time scale shorter than the period of oscillation. This is the regime in 
which viscosity erases all perturbations in the fluid. 

All these expectations can be checked in flgure 3 (right). For the reference model, the mode 
Ss{k = lO^^Mpc) experiences a very short stage of gravitational ampliflcation and is then driven to 
zero. The mode Ss{k = 10~^Mpc) has three visible stages: ampliflcation, oscillation (for just half- 
a-period) and damping. The mode (5s (fc = 10~^Mpc) is already in the acoustic oscillation phase at 
the earliest time shown in the flgure. Around r = 3000 Mpc, it also experiences a brutal viscosity 
damping. 

The conclusion of this section is that the viscous fluid perturbations never grow signiflcantly 
during matter domination and are washed out when viscosity becomes important (i.e. prior to dark 
energy perturbation). From this point of view, our viscous fluid model should be indistinguishable 
from any other dark energy model with negligible perturbations, like e.g. a scalar fleld with Cs = 1. 

The only time at which the viscous fluid perturbation can play a role is in the radiation domina- 
tion epoch. During that stage, the fluid dilutes like ps ex a^* and from the background point of view 
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X (Mpc) 



Figure 3. Left: Evolution of a few characteristic wavenumbers: that of Hubble crossing, kn = ciH (bottom 
black solid line); the Jeans wavenumber (upper blue solid line); and an effective Jeans wavenumber including 
some effect from viscosity (dashed line; see the text for details). Right: Evolution of the density perturbation 
5s for three wavenumbers between conformal time r = 200 Mpc and today. These two plots are obtained for 
the reference model with (e, A) = (0.7,0.25). 

it can be treated as extra massless neutrino species, as explained in section 3.2. This is not obvious at 
the perturbation level, since decoupled massless neutrinos are collisionless with anisotropic pressure 
leading to shear viscosity damping. Instead, our fluid is self-coupled and we assume its anisotropic 
stress to vanish. Hence, during radiation domination and before the time at which viscosity becomes 
important, we expect Ss to oscillate with a larger amplitude than (5,y, even if these two fluid share the 
same initial conditions and background evolution. This effect could in principle be significant. For 
instance, it has been proved that the damping of neutrino oscillations caused by shear viscosity is 
measurable with WMAP data [49] . If detectable, this effect would lead to slightly higher CMB peaks 
in the viscous model than in a wCDM model with the same A^cff- 

These expectations can be checked by running CLASS for a pair of models: a viscous fluid model 
and an extended wCDM model with the same values of (iVeff, wq, Wa) obtained with the fitting 
formulas 3.13-3.15. For such pairs of models and three choices of parameters (same as in section 3.2), 
we present in figure 4 the CMB temperature and matter power spectra. For the model with e = 0.1 
(A'off — 3.04), the matching is impressive and confirms our expectations. For the other two models 
with e ~ 0.7 (A'^^ff ~ 3.57), slightly larger CMB peaks can indeed be observed in the viscous model, 
for the reason explained above. This effect is however small, and we recall that e = 0.7 is an extreme 
assumption: it implies that the scalar particles decouple from the thermal bath late enough so that 
no entropy creation occured afterwards, except at the time of electron-positron annihilation. 

5 Compatibility with observations 

If the viscous fluid or "Dark Goo" model was leading to very specific predictions, we would perform 
a separate analysis to check its compatibility with observations and would infer confidence limits on 
each free parameters. However, we establish in sections 3.3 and 4.2 that at the level of background 
and linear perturbation evolution, the minimal Dark Goo model is nearly equivalent to an extended 
wCDM scenario with the following free parameters beyond those of ACDM: 



• An effective neutrino number 



A^cff = 3.04 + 2.2e^ 



(5.1) 



• A current value of the equation-of-state parameter 

Wo = -0.9085 + 0.21(r2„ - 0.3) + 3A"°-^10-* 



(5.2) 
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Figure 4. CMB temperature spectrum (left) and matter power spectrum (right) for pairs of models. In each 
plot, we present one viscous fluid model and a nearly equivalent extended wCDM model with the same values 
of {Ncff, Wo, Wa) obtained with the fitting formulas 3.13-3.15. From top to bottom, we assume one of the 
three models already studied in section 3.2. 



A time variation of the equation-of-state parameter 



0.129 -f 0.02{n,n - 0.3) + \-"nO 



(5.3) 



Various analysis have been carried for models with free N^e or (wo^Wa), but these parameters are 
usually not explored altogether at the same time. However, in our case the three quantities {Nes, 
wq, Wa) are expected to remain close to their minimal ACDM value (3.04, -1, 0), and by considering 
bounds obtained separately on these parameters, we do not expect to make a large error. We also 
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Figure 5. Predictions of the Dark Goo model for the parameters wo and Wa, compared to the 68.3%, 
95.4% and 99.7% confidence regions obtained from a joint fit of Supernovae, BAO and CMB data [52]. Left: 
Prediction for wq for each given Sim, when the other parameters of the model are varied in the following 
ranges: h £ [0.714,0.724], and A G [2 x 10"^, 1] (corresponding to mo G [10"^, 2] eV). Right: Predictions for 
Wo and Wa when other parameters are varied in the ranges described above and flm G [0.245, 0.306]. 



neglect as a first approximation the fact that the viscous model predicts slightly enhanced CMB peaks 
when Ncff is significantly increased (i.e. when e > 0.5). 

The strongest observational limits on A^eff come from Big Bang Nucleosynthesis, Ncff = 3.8 ±0.8 
at the 95% C.L. [50], while a combination of CMB and LSS experiments shows a marginal preference 
for values N^g > 3 [3, 51]. If this preference is confirmed, our model could account for any value of iVoff 
in the range [3.04, 3.62], as long as we assume that the scalar particles decouple before electron-positron 
annihilation such that e < T^/T^ (or eventually for N^s in the range [3.62, 5.25] if they decouple later 
such that e < 1, but this case would require a specific study of Big Bang Nucleosynthesis and photon 
recombination, and is not covered by the present paper). However, current error bars on NcS are large 
and do not allow to discriminate between models with a late scalar decoupling (3.04 < A'cff < 3.62) 
and early decoupling (A^cff — 3.04). 

Neglecting the effect of a possibly enhanced value of A^oft, we can readily use the bounds derived 
in the literature on the parameters {^m, wo, Wa) in order to test the compatibility of our model 
with current data. We first consider the region in (J7,„, wq) space preferred by a conrbination of 
Supernovae, BAO and CMB data presented in ref. [52]. If we assume h G [0.714,0.724] (as indicated 
by WMAP+BAO) and A G [2 x 10"^ 1] (the range of validity of our model), we find that (f^™, wq) 
must lie inside the thin blue band shown in the left panel of figure 5. Note that this band is obtained 
numerically with CLASS, but using our analytical fitting formula 3.14 would be equivalent. This band 
has a large overlap with the 68% preferred region of ref. [52] , thus the usual bounds on ftrn also apply 
in our case. This comparison does not allow to discriminate between different values of A, since the 
impact of varying A on wq is very small compared to the sensitivity of the data to wq ■ 

Finally, we vary {h, A) in the same ranges with r2,„ in his preferred range [0.245,0.306] and 
obtain an allowed region in {wq, Wa) space. This region is very small and reduces almost to a point 
(—0.9, 0.13) with respect to the region preferred by the data, as can be checked on the right panel of 
figure 5. This point lies in the 95% preferred region, very close to the 68% region, showing that our 
model is still compatible with current data. 

In the future, with better bounds on Ncs, wq and Wa, it will be possible to better test and nraybe 
to exclude our mininral Dark Goo model, especially if the values {wq, Wa) ^ (^0.9, 0.13) can be ruled 
out at some point. 
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6 Conclusion 

In summary we present a model of cosmological evolution with an additional scalar field. We show 
that in a certain range of values for the mass and self-coupling of the scalar field, bulk viscosity 
plays an important role at late times and mimics a dark energy behavior. At the background level, 
the model is compatible with current data on the equation-of-state parameter and predicts a small 
time variation for this parameter. At the perturbation level, the bulk viscous model produces the 
same temperature anisotropics and matter power spectrum as an extended wCDM model, except in 
the extreme case where the scalar field decouples from the photons at a very late stage and small 
deviations in the CMB peak heigths are observed. The model can also account for the extra rclativistic 
degrees of freedom that are marginally preferred by the data. 

The Dark Goo model has several features worth mentioning. For instance, all parameter ranges 
are bounded, making the model easily falsifiable by future dark energy experiments (note that the 
upper bound comes from the use of perturbation theory and is not a strict bound in that sense) . The 
possible values for the model parameters are "reasonable" and there is no apparent fine tuning. On 
a more philosophical level, it is also reassuring to know that the model has a built-in mechanism (i.e. 
breakdown of hydrodynamics) that may prevent the universe from accelerating forever (no "Big Rip" ). 
The model is also microscopic and the functional form for the bulk viscosity is obtained from first 
principles. This allows for a discussion of the validity of the model in terms of microscopic quantities 
(coupling, mean free path, etc). This is to be contrasted with other bulk viscous or Chaplygin gas 
models, where an exotic equation of state with little or no physical justification is postulated. 

We stress that the whole model hinges on the validity of hydrodynamics. The assumptions of 
hydrodynamics have to be satisfied in order for the concept of bulk viscosity to make sense. We 
discuss these issues and provide plausibility arguments for each assumptions, but these arguments 
can certainly be improved. In particular, a more fundamental understanding of cavitation would be 
helpful. 
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A Modifications to perturbation equations in CLASS 

In section 4, we present the perturbation equations for the viscous fluid using the conventions of 
ref. [47]. In order to make our results easily reproducible, we would like to recapitulate in this appendix 
the list of changes implemented in the CLASS code in order to account for these perturbations. CLASS 
uses the notations of ref. [48]. Density perturbations are represented by (J^'s and velocity divergences by 
0,;'s. We implement the new equations in the usual synchronous gauge comoving with cold dark matter 
{Ocdm = 0), with metric perturbations ij and h. The code uses conformal time r, and throughout this 
appendix the dots denote derivatives with respect to r, unlike in section 4. We define H = d/a = aH. 
When the system of perturbed equations is integrated over time for each wavenumber, at each new 
step we interpolate in the table of background quantities in order to obtain ps, Ps, Ci Pcff = Ps — SH(^, 
D = p'g, and finally 



dx 27r^ \ x'^ x^ J \ mo x 

We then infer the quantities Wg = Pcs/Ps and 

Ws = m{l + Ws){-Pcs/D + Ws) - 3HC/ps, (A.2) 
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(note that there is no typo here: this formula mixes % and H). The effective pressure perturbation 
5peS is derived from 

5p.« = ^PsSs - a-' U + ^ J C, (A.3) 

which is equivalent to cq. 4.11, replacing Sx by p^Sg/D, and with the identification ip = (dh/dt)/2 = 
a^^h/2 and a^^W^SU = Os- The continuity and Euler equations for the viscous scalar fluid read: 

Ss = -(1 + ws) (o, + ^ J - 3H (^ - WsS^ , (A.4) 

l+Ws {l+Ws)Ps 

In the Einstein equations, we include psSs and (ps + Pcff)Os in the calculation of the total Sp and 
ip + P)^ source terms. 

Finally, as far as initial conditions arc concerned, we treat the viscous scalar fluid at very early 
times as any other extra relativistic degrees of freedom. There is actually a difference between masslcss 
neutrinos and our fluid at early times: the latter are coUisionless, while the former are self-coupled. 
This implies that neutrinos can have a signiflcant anisotropic stress, while we assume that the scalar 
fluid can maintain an isotropic pressure. However, on super-Hubble scales on which initial conditions 
are implemented, this difference should have a negligible impact. 

Usually, the expression of initial conditions for all species involves a parameter R^, defined as the 
ratio of the masslcss neutrino density over the total ultra- relativistic matter densities R^ ~ Pv/ Pr- We 
include ps in both the numerator and denominator. We then derive initial conditions for all species in 
the usual way, including the density and velocity perturbations of ultra-relativistic species ((5ur, ^ur)- 
We initialize [Sg, 6s) to precisely the same values as [Sun ^ur)- 
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